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Abstract Pulsar timing is a technique that uses the highly stable spin periods of 
neutron stars to investigate a wide range of topics in physics and astrophysics. Pul- 
sar timing arrays (PTAs) use sets of extremely well-timed pulsars as a Galaxy-scale 
detector with arms extending between Earth and each pulsar in the array. These chal- 
lenging experiments look for correlated deviations in the pulsars’ timing that are 
caused by low-frequency gravitational waves (GWs) traversing our Galaxy. PTAs 
are particularly sensitive to GWs at nanohertz frequencies, which makes them com- 
plementary to other space- and ground-based detectors. In this chapter, we will de- 
scribe the methodology behind pulsar timing; provide an overview of the potential 
uses of PTAs; and summarise where current PTA-based detection efforts stand. Most 
predictions expect PTAs to successfully detect a cosmological background of GWs 
emitted by supermassive black-hole binaries and also potentially detect continuous- 
wave emission from binary supermassive black holes, within the next several years. 


J. P. W. Verbiest 

Fakultät für Physik, Universität Bielefeld, Postfach 100131, 33501 Bielefeld, Germany and 
Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, 53121 Bonn, Germany, e-mail: 
verbiest@physik.uni-bielefeld.ds 


S. Ostowski 
Centre for Astrophysics and Supercomputing, Swinburne University of Technology, PO Box 218, 
Hawthorn, VIC 3122, Australia, e-mail: stefan.oslowski@gmail.com 


S. Burke-Spolaor 

Department of Physics and Astronomy, West Virginia University, P.O. Box 6315, Morgantown, 
WV 26506, USA and 

Center for Gravitational Waves and Cosmology, West Virginia University, Chestnut Ridge Re- 
search Building, Morgantown, WV 26505, USA and 

Canadian Institute for Advanced Research, CIFAR Azrieli Global Scholar, MaRS Cen- 
tre West Tower, 661 University Ave. Suite 505, Toronto ON MSG 1M1, Canada, e-mail: 
sarahbspolaor@gmail.com 


* corresponding author 


2 J. P. W. Verbiest, S. Ostowski and S. Burke-Spolaor 


Keywords 


Pulsars; Pulsar Timing; Timing Array; Black Holes 


Introduction 


Neutron stars are the collapsed cores of massive stars that have undergone a su- 
pernova explosion after the end of nuclear burning and are supported from further 
collapse by neutron degeneracy pressure [12, 52, 114]. Since neutron stars are far 
more compact than their progenitor stars, they tend to exhibit very short rotational 
periods and extremely strong magnetic fields, as shown in Figure 1. Generally the 
magnetic axis is not aligned with the spin axis, so magnetic dipole radiation that is 
created in the neutron star’s atmosphere is swept around in space, somewhat like 
the beam of a lighthouse (this is the so-called “lighthouse model”). Depending on 
the orientation of the beam and its width, Earth may fall within that radiation beam 
once per rotation — which then causes the neutron star to be detected as a source of 
pulsed radiation, otherwise referred to as a “pulsar”. 


Radio Emission from Pulsars 


Following the lighthouse model described above, it is natural to expect that pulsars 
appear to the observer as so-called “pulse trains”: pulses of emission separated by a 
fixed period that equals the spin. These appear with a shape defined by the plasma 
properties in the pulsar’s magnetosphere, which can differ greatly from one pulsar 
to the next (see Figure 2). The emission mechanism of pulsars is understood in 
broad terms [see e.g., 109]. A few intriguing observational facts have been identified 
over the half century since the first pulsars were discovered. Most importantly, it 
has been shown that for most pulsars, the exact shape of individual pulses changes 
randomly from one period to the next. In contrast, however, the average shape of the 
pulsed emission is typically stable on timescales from minutes up to decades [56]. 
This implies that for any given pulsar, the pulsed emission can be averaged after 
accounting for the pulsar’s rotational period. The average pulse shape that can be 
obtained in this way is unique for the pulsar and can be thought of as a fingerprint. At 
radio wavelengths this average and reproducible pulse shape is typically called the 
pulse profile of the pulsar, whereas the term light curve is more common at higher 
frequencies (gamma and X-rays). The shape of the profile is defined by the emitting 
geometry in the pulsar magnetosphere. Since it is expected that the emission height 
is different for photons with different frequencies [36], it stands to reason that the 
shape of the pulse profile also typically differs with observing frequency (again, see 
Figure 2). 
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Fig. 1 P — P diagram for all pulsars known. Shown are the spin period and spin-period deriva- 
tive for all pulsars included in the ATNF pulsar catalogue version 1.61 [105]. Solitary pulsars are 
shown as full dots and pulsars in binary systems are shown as open circles. The grey dotted lines 
slanting downwards from the left represent surface magnetic-field strengths of 1017 G, 10!! G and 
10° G from top to bottom; and the black dot-dashed lines slanting upwards to the right represent 
characteristic ages of 10* yr, 108 yr and 10!° yr respectively, also from top to bottom. 


Not all pulsars have a stable pulse profile and not all pulsars emit radiation all the 
time. Indeed, a veritable zoo of pulsar emission phenomena has been discovered, 
studied and described throughout the years. There are so-called “nulling” pulsars 
[15], which often turn off, only to reappear at some point after. Some pulsars null for 
minutes on end, others for hours — some even turn off for months or years (at which 
point they are also called “intermittent” pulsars), suggesting an almost continuous 
distribution all the way up to so-called “RRATs” (Rotating Radio Transients), which 
only sporadically emit one or several pulses of radiation [108, 29]. Another category 
of pulsar emission is displayed by the “moding” pulsars [14]. These do not have one 
characteristic fingerprint, but two or three — and they arbitrarily change between 
them: while one day their pulse profile may look one way, the next day it may 
look different, only to go back to its original state on day three. Moding can also 
have a wide range of possible time scales, from single pulses all the way up to 
months or years between mode changes. Finally, there are drifting pulsars [46]. 
These “drifters” also have a well-defined pulse period that is readily and repetitively 
measurable on timescales of minutes to hours, but on timescales of seconds that 
pulse period seems to be overestimated, as the pulse appears to come a bit too late 
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Fig. 2 Pulse profile shapes vary across pulsars and observing frequencies. Shown are pulse pro- 
files for two traditional PTA pulsars: PSR J0437—4715 (left column), and PSR J1022+1001 (right 
column); and for each pulsar this profile is shown at three different observing frequencies: at 
~ 726 MHz (top row), ~ 1369 MHz (middle row) and ~ 3100 MHz (bottom row). Total intensity 
profiles are shown in full black lines, linear polarisation in dashed red lines and circular polari- 
sation in blue dotted lines. The tendency of getting sharper profile shapes at higher frequencies 
causes the timing precision to increase at those frequencies, although generally the noise level also 
increases due to the steep spectrum of pulsars [20, 65]. Consequently, most pulsar timing to date 
has been carried out at intermediate frequencies around 1.4 GHz. Finally, while this figure only 
shows profiles for MSPs, the evolution with frequency has been shown to be far more extreme in 
the case of slow pulsars [79]. (These plots were made based on the public data published by Dai 
et al. [39] and contain hundreds of hours of observing time at the higher frequencies, causing the 
radiometer noise to be barely visible.) 
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after each rotation, only to “reset” after a few dozen rotations, leading to the more 
stable long-term periodicity. 

Luckily for pulsar astronomers, the nulling, moding and drifting pulsars have 
turned out to be the exception rather than the rule and the majority of pulsars 
manifest themselves as predictably repetitive pulses of emission that have arbitrary 
pulse shapes on timescales of seconds, but stable and well-defined pulse shapes on 
timescales from minutes upwards. 


Pulsar Lifecycle and Spin Properties 


Whereas the first few pulsars that were discovered appeared to be a fairly homoge- 
neous group of objects, extensive surveys have continuously expanded the param- 
eter space in which pulsars have been discovered; and consequently a wide variety 
of pulsar types is now known. 

After the formative supernova explosion, pulsars start their new life as so-called 

with spin periods of a few tens of milliseconds and a magnetic field 
strength of order 10!” G at their surface [132, 24]. The emission of magnetic dipole 
radiation does make them lose angular momentum and consequently their rotation 
slows down gradually, typically by about 107}? seconds per second. The youngest 
and most well-known example of this class of pulsars is the Crab pulsar which was 
formed in a supernova explosion in 1054 CE [97, 136], which has been recorded by 
several civilisations across the World. 

After the first few thousand years, the spin period of young pulsars has slowed 
down sufficiently to have an appreciable impact on the spindown itself, which slows 
down their evolution. At this point they turn into the first discovered type of pulsar: 
the so-called “slow” pulsars or “canonical” pulsars. These are pulsars with rotation 
periods between about a tenth of a second and roughly ten seconds. Their mag- 
netic fields are thought to have strengths of roughly 10!°G to 10! G and they are 
expected to have formed thousands to hundreds of millions of years ago (see Fig- 
ure 1). Their spin period gets longer by about 107!4 s every second, as they lose 
angular momentum slowly but steadily. This loss in angular momentum causes them 
to eventually rotate too slowly to produce detectable amounts of radio waves and so 
after about a billion years they “turn off” and become undetectable [125]. 

Most of the slow pulsars are solitary objects that were either born as solitary stars, 
or broke free from their companion stars during their supernova explosion. A small 
sub-set, however, do have companion stars that are almost without exception main- 
sequence stars. When these main-sequence stars evolve and turn into red giants, it is 
not uncommon that their outer atmospheric layers stray into the gravitational well of 
the pulsar and cause it to accrete matter; and along with the matter, angular momen- 
tum. These pulsars are then spun-up while their magnetic field decays. The result is 
a “millisecond” pulsar (MSP), with spin periods between 1 ms and about 30 ms. Due 
to the accretion process, MSPs have relatively weak magnetic fields (10° G or less) 
and consequently their spin-down rates are also far lower than for slow pulsars (typ- 
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ically of the order of 107” s/s). As a result the rotational and emission properties of 
MSPs are not thought to appreciably evolve over their lifetimes. For an overview of 
formation and evolution of MSPs see [3, 22]. 

Given the extremely small spin periods, stable pulse profiles for MSPs can al- 
ready be obtained in a matter of minutes or even seconds. Furthermore, to date, only 
very few MSPs [101] have been demonstrated to show any of the anomalous emis- 
sion properties (nulling, moding, drifting) mentioned in the previous section that 
some of the slow pulsars display. Finally, due to the much larger angular momen- 
tum, MSPs have turned out to be far more stable clocks than slow pulsars. These are 
the reasons why MSPs have become known as “nature’s gift to physics”: the perfect 
Einstein clock that can be used to test a wide range of relativistic predictions. 


Pulsar Timing and Pulsar Timing Arrays 


Pulsar timing is a method that exploits the highly regular spin period of pulsars and 
their predictable pulse shapes, to study a wide range of questions in physics and 
astrophysics. In essence, when doing pulsar timing, one monitors the times at which 
subsequent pulses from a pulsar arrive at an observatory. These observed pulse-ar- 

or ToAs are then compared to a mathematical model that attempts to 
quantify all the factors that impact the travel time of the electromagnetic waves on 
their way from the pulsar to the Earth. In practice, a number of complications need 
to be dealt with, as outlined below. A more advanced approach is to use multiple 
pulsars — an “array” of pulsars — to look for signals that correlate between different 
pulsar pairs. Such experiments are called “pulsar timing arrays” (PTAs). 


Template Profiles 


In its simplest form, pulsar timing could be based on the times when the peaks in 
a train of pulses are detected. In order to increase the measurement precision, one 
could also take the intensity-weighted average arrival time of any given pulse. A far 
more powerful method, however, is to use the information encoded in the shape of 
the pulse, to measure the arrival time relative to a standardised pulse shape. This can 
be thought of as the ultimate, noise-free pulse profile. Obtaining such a standardised 
pulse shape or “template profile” is necessarily an iterative process. Fundamentally, 
as many pulses should be averaged together as possible. However, in order to align 
said pulses, an accurate pulsar timing model should be used to predict the phase of 
subsequent pulses to a precision far better than the time-resolution afforded by a 
typical observation. To give an example, we aim to predict arrival times to a preci- 
sion of nanoseconds while typical observations have phase resolution on the order 
of microseconds. The use of all of the information contained in a complex pulse 
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profile, through the use of a template profile, allows one to achieve this necessary 
magnitude of improvement. 

Since the creation of the template profile itself is the very start of the path towards 
a functional pulsar timing model, typically the entire process gets iterated so that 
the template profile and the timing model can both improve until their solutions 
converge. 

A danger in the creation of template profiles is so-called “self-standarding”’. This 
is a phenomenon that occurs when the data that are being used in the timing, are 
also used to construct the template profile. Specifically, if the number of observa- 
tions that are combined to construct the template profile is too low (a typical rule 
of thumb is that “too low” is less than a few thousand), then the noise within the 
observation can be “recognised” in the template profile, leading to inaccurate offset 
measurements [as illustrated clearly in Appendix A1 of 64]. Consequently, many 
timing experiments make use of analytic models to describe the template profile. 
These may not always be able to perfectly model all the pulse-shape features, but 
they avoid timing corruptions caused by self-standarding. Alternatively, a smooth- 
ing filter may be applied to the template, in order to reduce the correlating noise 
[41]. 

It was mentioned earlier that pulse profile shapes typically change with the ob- 
serving frequency. This should ideally be taken into account when constructing the 
template profile, i.e. the template profile should effectively have a dependence on 
observing frequency, too. In most published works this has not been the case be- 
cause the bandwidth of pulsar observations used to be sufficiently narrow that fre- 
quency evolution of the pulse profiles was effectively undetectable, but over the 
last decade (fractional) bandwidths of observing systems have increased so signif- 
icantly that so-called “frequency-resolved” template profiles are rapidly becoming 
the norm rather than the exception. Also in this case, one can either create an ana- 
lytic description of the profile in two dimensions [1 15, 95] or use a template purely 
based on accumulated data [45]. 


Template Matching 


Once a template profile has been created, it can be used to calculate the ToAs of the 
various observed pulse trains. Since most pulsars are so faint that individual pulses 
cannot be detected and because single pulses are usually not all alike, standard 
pulsar-timing experiments do not time individual pulses, but average subsequent 
pulses modulo the pulse period!. This averaging procedure is commonly referred to 
as folding — it reduces the time resolution of the observations while phase resolu- 
tion is maintained. Typical values for the time resolution of pulsar-timing data after 
folding, are anywhere from minutes to one hour, depending on the brightness of the 
pulsar and the goal of the experiment. Phase resolution is defined by the number of 


1 At this point the question of which pulse is being timed exactly, is in principle an arbitrary choice, 
but the most commonly used pulsar-timing software uses a pulse in the centre of the observation. 
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bins across the profile, with typical values ranging from 128 to 2048, depending on 
the sensitivity of the telescope and the bluntness or sharpness of the pulse profile. 

The folded pulse profiles could be cross-correlated with the template profile in 
order to achieve the phase of the observation — which can then be added to the ob- 
servation’s time stamp in order to achieve a ToA. In practice this measurement is 
commonly undertaken in the Fourier domain, as explained in detail in the appendix 
of Taylor [143]. Since an offset in the cross-correlation is equivalent to a phase gra- 
dient in the cross-power spectrum, typically ToAs are determined by least-squares 
fitting the phase gradient in the cross-power spectrum of the template profile and 
the folded observation. The phase offset resulting from this is then added to the 
observation’s time stamp. 

The measurement uncertainty of these phase offsets — and hence of the ToAs — is 
an important value as well, since many pulsars appear to have highly variable flux 
densities (a process caused by the interstellar medium, called scintillation), which 
means that not every ToA carries as much information and hence should not be 
weighted equally. Specifically, the pulse profile that is to be timed will contain a 
certain amount of white noise called radiometer noise, which depends on the obser- 
vational properties of the pulsar and the observing system as follows [96]: 


Oroa = S/N = Bue FE ELE a) 
sys 

where Oroa is the ToA uncertainty, S/N is the signal-to-noise ratio of the observa- 
tion, B is a factor which describes instrumental losses, e.g. due to digitisation, np, 
tin. and Af are respectively the number of polarisations combined, the integration 
time and the bandwidth of the observation, Tpeak is the brightness temperature of 
the pulsar at the peak of its profile and Te is the brightness temperature (i.e. noise) 
of the observing system, which typically includes corrections for cable losses, in- 
strumental gain, spill-over and sky noise, amongst other things. W is the equivalent 
width of the pulse profile, defined as the integrated pulse intensity divided by the 
peak intensity and P is the pulse period. 

The radiometer noise is the most fundamental factor limiting pulsar-timing pre- 
cision, in the sense that it is present in all observations and is determined to a large 
degree by the fixed properties of the pulsar and the technical capabilities of the tele- 
scope. Traditionally it was quantified as the formal uncertainty of the phase-gradient 
fit described above, although it has been shown that in the low-S/N regime this can 
cause irregularities [see 7, App. B], leading people to either remove ToAs below a 
certain S/N level (e.g. requiring S/N > 8) or to apply more advanced, Monte-Carlo- 
based uncertainty estimations, as was proposed as “good pulsar timing practice” by 
Verbiest et al. [153]. 
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Timing Model Determination 


Once the ToAs have been measured, they need to be compared to predicted arrival 
times provided by a pulsar timing model. A timing model is a mathematical formula 
that predicts the arrival time of a pulse based on a set of timing-model parameters. 
Generally, the timing model is defined in two steps [see also 143, 47]. In the first 
step, the measured pulse arrival time fops is referred to a time of emission at the 
pulsar tpsr. This is achieved by accounting for all known propagation and geometric 
delays: 
ÍPSR = fobs — Ao — Attsm — Apin- 


Specifically, first the pulse ToA is transferred to the Solar System barycentre (i.e. 
corrected by a delay Aj), which is the inertial reference frame most commonly used 
in pulsar timing. This transformation includes correction factors for relativistic ef- 
fects caused by the mass distribution in the Solar system and the Earth’s orbital and 
rotational velocity, for atmospheric propagation delays (note these have mostly been 
neglected to date but will become important with the next generation of radio tele- 
scopes), for light-propagation times (the so-called Roemer delay), parallax effects, 
frequency-dependent propagation delays induced by the Solar Wind and corrections 
for the observatory clock. 

After the ToAs have been transferred to the Solar System barycentre, interstel- 
lar propagation delays (Aysy) are corrected for. As described in more detail in the 
next section, these delays have long been treated as dependent on the observing fre- 
quency, but constant in time. With increased measurement precision provided by 
wider bandwidths and lower observing frequencies, time-variable models of inter- 
stellar dispersive delays are now becoming more common. 

For pulsars that inhabit binary systems, there is one further transformation, 
namely from the barycentre of the binary system to the pulsar (delay Agin). This 
includes the Roemer delay based on a Keplerian description of the binary orbit, but 
can also contain relativistic effects such as the Shapiro Delay (time dilation caused 
by the companion star’s gravitational field), the Einstein delay (time dilation caused 
by the pulsar’s gravitational field and gravitational redshift) or a host of other more 
complex effects, dependent on the binary’s properties. [See 47, for a complete list- 
ing.] 

Once the time of emission fpep is determined, it can be converted to a rotational 
phase based on a spin-down model that is usually simply described as a Taylor 
expansion: 


1 
$ (tpsr) = V (fesr — t0) zf (fesr = to)” +--+ 2) 


where v is the spin frequency of the pulsar, V its first derivative and tọ an arbi- 
trary reference epoch. Standard electromagnetic theory predicts a second frequency 
derivative V = av but in practice this is immeasurably small in the case of MSPs 
and is typically obscured by other effects (so-called timing noise, see further) in 
most slow pulsars. Consequently, by default pulsar timing models contain a spin 
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frequency and frequency derivative but not usually any higher-order spin frequency 
derivatives. 

Initial timing models are derived from pulsar-search observations. These are raw 
time series that are not folded, but instead are Fourier transformed in order to ob- 
tain an instantaneous pulse period. By monitoring the pulse frequency evolve over 
several such observations, an initial estimate of the spindown and of the binary orbit 
can be determined. This then constitutes an initial timing model that can be used to 
predict the pulse frequency for future observations, allowing the data to be folded in 
real-time, which makes observations much less demanding in terms of data-storage 
and processing power requirements. 

When a basic timing model has been constructed that is able to predict the arrival 
time of future observed pulses to well within a pulse period, it is said that “phase 

"has been achieved. From this point forward, the phases calculated in 
Equation 2 can be used to improve the timing model. Effectively, these phases 
should all be zero if the timing model was perfect and no corrupting noise sources 
were present. Consequently, any deviation from zero highlights effects which do af- 
fect the observations, but are not taken into account (correctly) by the timing model. 
These differences between the observation and the model are called the 
uals and they lie at the core of pulsar-timing analyses. Indeed, the art of pulsar tim- 
ing is to optimise and extend the timing model to decrease these timing residuals. 
In the optimal case, the timing residuals will only consist of white noise which is 
accurately quantified by the uncertainties of the ToAs. In this case the timing is fully 
dominated by the radiometer noise described earlier. 

Each parameter that is contained in the timing model affects the timing residuals 
in a well-defined way, which is called the timing signature of this parameter. Con- 
sequently, simply by visual inspection of the timing residuals, particular errors can 
sometimes easily be picked out, as shown in Figure 3. In most practical cases the 
uneven spacing between observations, the variability in the ToA uncertainty and the 
combination of multiple timing signatures in incomplete or outdated timing mod- 
els make the analysis of timing residuals rather more complex than these simplified 
examples suggest. 


Interstellar Propagation Delays 


A particular source of difficulty when analysing pulsar-timing data, is the impact 
of the ionised part of the interstellar medium (also referred to as the ISM). The 
refractive index of the interstellar medium is determined by the plasma frequency, 
which is dependent on the local electron density [96]: 
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Fig. 3 Examples of pulsar timing residuals. The top-left figure shows a typical PTA data set taken 
from [106]. The residuals are not purely white in this case, most likely due to a combination of 
uncorrected variations in interstellar propagation delays and interstellar dispersion. In the top-right 
figure the impact of a 1% change in the spindown is demonstrated, leading to a clear quadratic trend 
in the residuals as the spin period gets increasingly incorrect as time progresses. The bottom-left 
figure shows a positional offset of 0.1 arcsec in both right ascension and declination, leading to an 
annual sine wave with constant amplitude. The bottom-right plot shows what happens in contrast 
if the position is correct (at the reference epoch near the start of this data set) but the proper motion 
is 10% incorrect. This causes a linear increase in positional error and hence induces an annual sine 
wave with linearly growing amplitude. 


where e is the electron charge, ne is the electron density in units of cm~? and me 
is the electron mass. Given that plasma frequency, the refractive index for a photon 


2 
with frequency f can be determined as: u = 4/1 — (4) . Since the group veloc- 


ity of electromagnetic waves is dependent on the refractive index (vg = uco), this 
leads to a frequency-dependent group velocity which is observable as a frequency- 
dependent propagation delay: 


F- 


—2 
At = h DM, 
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where DM is the “dispersion measure” defined below and the constant K = 4 = 
2.41 x 10-4 MHz ?pc/em~?/s is the inverse of the “dispersion constant”. Theoreti- 
cally, the dispersion constant could be determined to higher precision (specifically, 
EE Sëch but in pulsar timing it has traditionally been fixed as given above [80]. 

The dispersion measure is straightforwardly defined as the integrated electron 
content along the line of sight between the telescope and the pulsar: 


D 
DM = f nedi, 
0) 


where ne is the electron density in cm~?, D is the pulsar distance in pe and DM is 
typically expressed in units of pc/cm7?. However, given that in most pulsar-timing 
software 2 has been defined fixed at the above-mentioned value and since the ac- 
tual observable is YxDM, in practice most DM measurements would require a 
slight correction before being interpreted as physical electron density measures, as 
described by Kulkarni [80]. 

Due to the typically high spatial velocities of radio pulsars [up to 1000 km/s and 
beyond 33], the line of sight along which the pulses travel to Earth sweeps through 
interstellar space; and due to the numerous turbulent structures that are present 
throughout that space [5], this motion causes the DM to change as a function of 
time. While such variability has long been known to exist, it was mostly detectable 
as a slowly-varying, red-noise process. Since the turn of the century, however, with 
improved instrumental sensitivity, ever more precise measurements of DM varia- 
tions in time have been detected and accurate modelling of DM(t) is becoming a 
complex and essential part of pulsar timing experiments [70, 74]. 

One particular component that contributes to temporal variations in DM, is the 
Solar Wind. Due to the annual motion of lines of sight — particularly for pulsars 
near the ecliptic plane — the additional dispersion caused by the Solar Wind has a 
clear annual signature, which is typically modelled straightforwardly by assuming 
the Solar Wind to be homogeneous and spherically symmetric [47]. Generally, it is 
assumed that such a straightforward model would suffice for the purposes of high- 
precision pulsar timing, except perhaps closest to the Sun, so in addition to the 
spherical models, PTAs have tended to remove ToAs for observations that took place 
within 5-10 degrees of the Sun [153]. It has been attempted to extrapolate optical 
observations of the Sun to derive a more detailed, inhomogeneous model of the Solar 
Wind for pulsar-timing purposes [157], but while this model was shown to provide 
an accurate spectrum of inhomogeneities, it does not provide accurate corrections 
for pulsar-timing experiments [146]. 

In addition to dispersion, the IISM introduces a host of other propagation effects, 
as recently reviewed by Stinebring [138]. In most cases these effects are not limiting 
timing precision yet, although time-variable scattering (also referred to as multi-path 
propagation — a phenomenon that widens the pulse shape through increased travel 
path lengths) has been shown to be relevant in the timing of at least one MSP [92]. 
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Timing Noise 


Probably the hardest effect to mitigate in pulsar timing is the so-called “timing 
noise”. This term is generically applied to any timing residuals that are not white 
noise and cannot be corrected for by any of the deterministic timing-model param- 
eters or by frequency-dependent DM models. Presumably this typically long-term 
noise is caused by inherent rotational instabilities in the neutron star itself [78], 
although the physical mechanism is as yet not known. 

Timing noise has been studied extensively in slow pulsars [99, 59], where it is 
highly common. In MSPs, timing noise has been shown to be far less common, or 
to exist only at much lower levels [152, 91]. Nevertheless, as the length of pulsar- 
timing data sets grows and the timing precision increases, the prevalence of timing 
noise — and the importance of mitigation techniques — continues to increase even in 
MSP timing projects [10]. 


Other Noise Sources 


After the timing model is optimised and the IISM variations are modelled and cor- 
rected for, ideally the timing residuals should be spectrally white and normally dis- 
tributed. In practice a wide range of effects can negatively affect the timing, as re- 
cently reviewed by Verbiest and Shaifullah [151], although in practice the primary 
impact aside from the IISM and timing noise, is pulse phase jitter, also known as 

. The work by Lam et al. [81] on 37 MSPs, shows that jitter is relevant pri- 
marily at low frequencies (particularly at observing frequencies below 1 GHz) but 
less so at higher frequencies, where high-precision pulsar timing is most commonly 
done. Since the importance of pulse jitter is strongly dependent on the sensitivity 
of the telescope (and on the available bandwidth), jitter will become an even more 
important source of noise in the next generation of radio telescopes [94]. One ap- 
proach to avoid jitter-dominated timing, would be to divide up large interferomet- 
ric telescopes into less-sensitive sub-arrays which allows more effective scheduling 
with lower jitter noise [135]. Another approach is to mitigate jitter noise in post 
processing as demonstrated in Osłowski et al. [112, 113]. 


Pulsar Timing Software 


All of the analysis described above tends to be carried out with two separate types of 
data-analysis tools. Firstly one needs software that can process the raw observation 
files: the so-called archives that store radio-wave intensity as a function of polari- 
sation, pulse phase, frequency and time. Secondly one needs model-fitting software 
that analyses the ToAs and the related timing models. 
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The primary software package that is used globally to analyse pulsar archives 
in the context of pulsar timing, is PSRCHIVE [63, 150]. The only exception to the 
use of PSRCHIVE for the analysis of pulsar archives and the creation of ToAs is the 
possible creation of broadband ToAs with analytic, frequency-dependent templates. 
For this purpose, the purpose-built PULSEPORTRAITURE software [115]? is used in- 
creasingly commonly. For the analysis of ToAs and timing models, however, a larger 
variety of software packages has been developed. The most common ToA-analysis 
package is TEMPO2 [62]*, which is fundamentally a C/C++ translation of the much 
older, Fortran 77-based TEMPO package [1 11]°, which is still being used, often in 
parallel with TEMPO2. More recently, the PINT package was developed [98]°, which 
is an independent timing package written in Python and which is primarily used in 
North America. 

Extending pulsar-timing software to constrain or detect correlated signals, such 
as those from GWs, is a non-trivial effort. Whereas frequentist methods have oc- 
casionally been implemented as part of packages like TEMPO2, it has become far 
more common to build independent software packages that are specifically aimed at 
Bayesian analyses of the timing model — including correlated signals. These pack- 
ages tend to be written in Python and use Python wrappers around the source code of 
the standard pulsar-timing packages mentioned above — most commonly TEMPO2. 
The most recently used package for such advanced ToA analysis (including GW 
analyses) is ENTERPRISE [48]’; two commonly used earlier packages are TEMPON- 
EST [89] and PICCARD [149]. 


Gravitational Waves and Other Correlated Signals 


In the previous sections we focussed on pulsar-timing effects and phenomena that 
affect each pulsar independently in fully unrelated ways. This is typical of the tra- 
ditional single-pulsar timing experiments that have been the mainstream of pulsar- 
timing research ever since the first pulsar discovery. However, thoughout the 1980s 
the realisation developed that some signals might have timing signatures that corre- 
late between pairs of pulsars [57, 123, 50]. These signals would require a new, more 
complex analysis as they require a joint analysis of a larger number, say an array, of 
pulsars. The concept of the PTA was born, but would not come to full fruition until 
the start of the new millennium, after the number of known pulsars (and the num- 


? http: //psrchive.sourceforge.net 

3 https: //github.com/pennucci/PulsePortraiture 
4https://bitbucket.org/psrsoft/tempo2/ 
Shttps://github.com/nanograv/tempo 
®https://github.com/nanograv/PINT 
Thttps://github.com/nanograv/enterprise 

8 https: //github.com/LindleyLentati/TempoNest 
H https://github.com/vhaasteren/piccard 
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ber of PTA-worthy pulsars) had dramatically increased following a couple of highly 
successful surveys, primarily in the inner Galaxy [103, 104]. In the following para- 
graphs, an overview will be given of typical correlated signals in pulsar-timing data, 
with particular focus on the signature of gravitational waves. 


Correlated Signals in Pulsar Timing Data 


As described above, pulsar timing is effectively a model-fitting exercise where deter- 
ministic parameters get optimised as increasing amounts of data constrain the timing 
model. In addition to the deterministic components of the timing model, there are 
non-deterministic effects like DM variations or jitter noise that affect the timing 
differently for each pulsar. Finally, there are signals — both deterministic and non- 
deterministic — that affect all pulsars in similar ways. Three types of such correlated 
signals have been described to date [123, 50]: 


e a monopolar signal that would most likely arise from imperfections in the refer- 
ence clock [60, e.g.], 

e a dipolar signal that is most likely related to the Solar-System ephemerides [53, 
e.g.] and 

e a quadrupolar signal that is predicted to arise from gravitational waves [57]. 


These correlated signals require more complex approaches: Solar-System ephemeris 
models could be updated with the pulsar-timing data [32, 53, 147], but clock signals 
and gravitational-wave signatures are in essence random and hence the actual cor- 
relations between timing of different pulsars need to be used in order to determine 
those. For monopolar (or clock) signals, this has been successfully achieved a num- 
ber of times [122, 60], but quadrupolar (or gravitational wave) signal have so far not 
been unambiguously detected, partly also because all these signals interact, making 
a clear detection of the highest-order correlated signal (i.e. the gravitational waves) 
dependent on accurate determination of all other correlated signals [145]. 


Effect of Gravitational Waves 


The effect of GWs on pulsar timing was first described by Detweiler [44] and more 
recently clearly summarised by Sesana and Vecchio [130]. Fundamentally, a GW 
passing over the Earth-pulsar system will introduce a time-variable redshift into the 
pulsed signal: 

z(t,Q) = (v(t, Â) = vo) Jun, 


where Q is the direction of propagation of the GW, Vo is a reference frequency 
and v(t, Q) is the observed pulse frequency which is defined by the geometry of the 
system (pulsar and GW position with respect to the observer) and the GW properties 
(polarisation and amplitude). The integral of these redshifts quantifies the impact on 
the timing residual: 
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for an observation taken a time t after the first observation in our data set. Further- 
more, the observed redshift is only dependent on the perturbation of the space-time 
metric at the position of the pulsar at the time when the pulse was emitted; and the 
space-time perturbation at the location of Earth when the pulse was received. 

A more straightforward way of putting this, is that the GW impact on the pulsar 
timing residuals has two equally large components: a so-called “pulsar term” that 
quantifies the impact of the GW on the emission of the pulse; and an" 
which quantifies the impact of the GW on the detection of the pulse on Earth. For 
non-evolving, sinusoidal GWs, the Earth and pulsar terms will be identical except 
for a phase offset between the two. Furthermore, the Earth term will affect all pul- 
sars equally (albeit modulated in a quadrupolar way, as described by the so-called 
Hellings-and-Downs curve [57], see Figure 4) whereas the pulsar term will have a 
different phase offset for each pulsar, since the phase of the pulsar term depends on 
the distance to the pulsar. 

For such mono-chromatic signals, this phase offset could in principle be mea- 
sured along with the pulsar distance, allowing extremely precise localisation of 
both the GW source and the pulsars in the array [87]. Such an experiment would, 
however, require a level of timing precision that is not realistically achievable with 
present-day telescopes (but may be achievable with the next generation Square Kilo- 
metre Array or SKA). At present, therefore, the pulsar term is typically considered 
a noise term, while the Earth term is the actual quadrupolar signal we hope to de- 
tect. The amplitude of the correlated signal will give insights into the GW’s origins 
(see the next Section), whereas the shape of the correlation curve could constrain 
fundamental physics of gravitational waves, like their polarisation properties [86] 
and propagation speed [84]. Furthermore, deviations from the theoretically expected 
Hellings-and-Downs curve can be expected for anisotropic backgrounds of GWs, 
where a few bright sources stand out above the background and cause a correlation 
function that is not only dependent on the angle between pulsars, but also on their 
location on the sky [144]. 


GW Sources in the PTA Band 


The sensitivity of PTAs to GWs is limited in terms of the GW frequency by the 
length of the observing time span and the cadence of the observations. Specifically, 
PTAs are most sensitive near a frequency of 1/T where T is the length of the data 
set, i.e. on the order of a decade or more, which corresponds to a frequency on the 
order of nanohertz. Since the GW impact is a change in the pulse frequency [44], 
but the timing residuals are phases, i.e. effectively integrals over pulse frequency, 
the sensitivity of PTAs decreases with increasing frequency down to the Nyquist 
frequency 2/C where C is the cadence of the observations, typically of the order of 
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Fig. 4 Correlated impact of a stochastic background of GWs on pulsar-timing residuals. Shown is 
the so-called Hellings-and-Downs curve [after 57] as solid black line, which quantifies the corre- 
lation induced into pulsar-timing residuals of sets a pulsars as a function of the angular separation 
between those pulsars on the sky. For a single GW the shape would be similar, but would be de- 
pendent on the orientation of the pulsar pair with respect to the GW source; for a stochastic back- 
ground only the angle between the pulsars matters. Square points show simulated measurements 
for an optimistic realisation of a PTA experiment. 


weeks or a month, which corresponds to a frequency cut-off of the order of micro- 
hertz. 

This frequency range makes PTAs particularly complementary with other GW 
detectors like LIGO [10 Hz-10kHz, 107] and LISA [0.1 mHz-1 Hz, 4]; and com- 
petitive with proposed GW detection methods based on space-based VLBI [28, 26]. 
The difference in GW frequencies furthermore implies that different sources of GWs 
can be expected to be detectable with PTAs. Specifically, four types of sources are 
anticipated, as described below. Example timing residuals induced by these four 
types of GW are shown in Figure 5. 


Gravitational wave background (GWB): A GWB is a superposition of GWs from 
a large number of GW sources that add incoherently. The most likely GWB in the 
PTA frequency range is widely expected to arise from 
binaries and predictions for its spectrum and amplitude are based on simulations 
such as the Millennium Run [129] or the Illustris simulation [75]. Specifically, 
the power spectrum of this GWB is expected to have a power-law shape with 
slope —2/3 and is likely to flatten or even tip over at GW frequencies lower 
than ~ 0.1 yr! [34]. Alternative backgrounds have been proposed [see, e.g. 27, 
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Fig. 5 Example timing residuals for four GW types on three different pulsars. These four panels 
show what the timing residuals that are caused by GWs could look like. The top-left panel shows 
the influence of a GWB (with characteristic strain amplitude 10~!> and spectral index —2/3), 
the top-right panel that of a CW (from a 10? Mo equal-mass binary supermassive-black hole at 
redshift z = 0.01), the bottom-left panel shows a BWM (with strain amplitude 5 x 10715) and 
the bottom-right panel shows a GW burst (without memory and with arbitrary waveform). The 
three different colours show the impact on different pulsars (i.e. different sky locations), with red 
showing the timing residuals of PSR J0437—4715, blue those of PSR J1012+5307 and black those 
of PSR J1713+0747. The simulated measurement uncertainty is 1 ns and no intrinsic spin noise 
or DM variations were included in the simulations. The top figure in each panel shows the pre-fit 
residuals (i.e. the raw GW signature), whereas the bottom plot shows the post-fit residuals. The 
difference between these is caused by fitting of the standard timing-model parameters. As can 
be seen, most of the power of the GWB is absorbed in the timing-model fit, since its long-term 
signature resembles the timing signatures of the pulse period and spindown, which absorb most of 
the signal. Such absorption of GW residuals in common timing-model parameters is less likely to 
happen for CWs or bursts, unless they happen to have a periodicity that is close to a year (or an 
integer fraction or multiple thereof) or to the orbital period of the pulsar being timed. 
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128, 127], with differences in both spectral shapes and amplitudes, implying that 
an eventual detection would be able to differentiate between the origins of the 
GWB, or would be able to place constraints on the galaxy-merger history of the 
Universe. 
Continuous waves (CWs): Continuous GWs in the PTA band are expected from 
single supermassive black-hole binaries that are close enough to Earth to stand 
out beyond the GWB. With improved sensitivity, a detection of CWs is generally 
expected to follow within a few years from a GWB detection [124, 110]. 
GW bursts: Single GW events like close encounters between supermassive black 
holes or some interactions [40], could result in a single burst of 
GWs, which might be detectable provided the burst itself lasts sufficiently long 
for it to affect multiple subsequent pulsar observations (i.e. at the very least days 
long); and provided the burst is sufficiently bright to stand out of the noise. 
Bursts with memory (BWMs): Bursts that are too faint to be detected directly, 
could still be detected as a “memory event” or a BWM. In this case it is the 
permanent deformation of the space-time metric [49] that has a lasting impact on 
the pulse frequency, causing the impact on the timing to accumulate over time. 


Present PTA Constraints 


Around the world, collaborations have emerged to carry out the large-scale observa- 
tional campaigns required for GW detection with PTAs. The Australian Parkes PTA 
[PPTA, see 106] was the first one, commencing observations in 2005 and revising 
some of the original theoretical work [67, 68]. Centred around the 64-m Parkes radio 
telescope, they have been monitoring between 20 and 30 pulsars using both dedi- 
cated observations [76] and archival data [152]. The PPTA last placed a limit on the 
GWB in 2015 [131], finding that the normalised amplitude at a frequency of 1 yr! 
must be lower than 10~!> with 95% confidence — which started to get into the area of 
theoretically predicted amplitudes at the time; this limit is still the most constraining 
bound on a GWB in the PTA band to date. Most recently, the PPTA has focussed 
its efforts on alternative sources of GWs, such as ultra-light scalar-field dark matter 
[118] as well as CWs and BWMs [100]; alongside a more intensive commitment to 
the global International PTA (see below) and instrumental development [61]. 

The European PTA [EPTA, see 43] also commenced in 2005, soon after the PPTA 
and also relies on a combination of specific PTA data and archival monitoring data, 
adding up to a total of 42 MSPs [43]. To date, the EPTA has primarily used data from 
the Jodrell Bank 76.2-m Lovell telescope, the 100-m Effelsberg Radio Telescope, 
the 300-by-35-m decimetric radio telescope at Nançay observatory and the Wester- 
bork Synthesis Radio Telescope which is an interferometric array consisting of 14 
antennae of 25-m diameter. Their most recent limit also dates back to 2015 [90], but 
is slightly less constraining, at Aj /y, < 3.0 x 107 15. The same data set has been used 
to place constraints on CWs [13] and possible anisotropies in the GWB [144]. Fi- 
nally, dedicated, high-cadence, data were used to place constraints in the microhertz 
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regime [116]. In recent years, the EPTA has primarily focussed on instrumental de- 
velopment with new data-recording systems [83, e.g.], the commissioning of a new 
64-m radio telescope in Italy [SRT, see 119] and the interferometric combination 
of all mayor EPTA telescopes in the “Large European Array for Pulsars” (LEAP) 
project [19]. The EPTA has also had a strong involvement in the commissioning 
of the LOw-Frequency ARray [LOFAR, see 148], which has shown to be a useful 
telescope for monitoring not only MSPs [77] but particularly DMs [45, e.g.]. 

The third major PTA is the North-American Nanohertz Observatory for Gravi- 
tational waves [NANOGrav 41]. NANOGrav has been particularly involved in pul- 
sar searches and as such their source list has been continuously expanding, count- 
ing 48 MSPs at their current data release [https://data.nanograv.org 
and 82]. The latest and most constraining limit on the amplitude of the GWB, is 
Ajjyr < 1.45 x 107" [9], slightly above the PPTA limit. Within the last few years, 
NANOGrav has also placed significant bounds on BWMs [2] and on CWs [1]; and 
placed a specific limit on the proposed binary black-hole system in the radio galaxy 
3C66B [142, 11]. 

The three original PTAs mentioned above joined forces to further increase sen- 
sitivity and in an attempt to decrease the time to the first GW detection in the nHz 
regime. As described by Manchester and IPTA [102], the first joint PTA meeting 
took place in 2008, but the formal establishment of the International PTA (IPTA) 
did not happen until 2011. Since then, the IPTA has released two combined data 
sets [153, 117] and has indicated that an improvement in GWB sensitivity by a 
factor of about two should be expected, but no full GW analysis has been carried 
out on IPTA data so far, given the complexities involved with the highly inhomo- 
geneous nature of the data set. These inhomogeneities and complexities [discussed 
and listed in 153] are a specific challenge for any combined PTA experiment and 
often requires additional research regarding detailed aspects of the analysis, or de- 
velopment of new methods that are able to deal with such inhomogeneous data. A 
lot of progress has been made in recent years in this regard, particularly with the 
advent of Bayesian analysis software packages like TEMPONEST [89, 91] and EN- 
TERPRISE [48], amongst others. So far, the first IPTA data combination has been 
used by Hobbs et al. [60] to construct a pulsar-based time standard (thereby solving 
for any monopolar correlations in the data), while Caballero et al. [31] used it to 
constrain errors in the Solar-System ephemeris models used. A comprehensive GW 
analysis is planned for the second data release. Meanwhile, analysis tools are being 
tested on mock data challenges [54, 151, 17]. 

As a new generation of telescopes is being constructed on the pathway to the 
SKA, a number of new PTAs have recently emerged. Specifically, the Indian PTA 
[InPTA, 71] has been formed in 2018 and uses high-precision timing data from 
the upgraded Giant Metre-wave Radio Telescope (aGMRT) and low-frequency data 
from the Ooty Radio Telescope (ORT). The Chinese PTA [CPTA 85] uses the Five- 
hundred meter Aperture Spherical Telescope (FAST), the Xingjiang Qitai 110-m 
Radio Telescope (QTT) and a network of 100-m-class radio telescopes across China 
and predict to go well below current sensitivity limits after even a few years of 
observing. Finally, the South-African MeerKAT telescope [69] is being used by 
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the international MeerTIME consortium [16] to produce (amongst other things) a 
highly sensitive PTA data set, part of which will be taken at relatively high radio 
frequencies, between 1.7 and 3.5GHz, thereby limiting the impact of interstellar 
effects. 


Table 1 Summary of present limits on GWs from PTA data. 


GW Type EPTA NANOGrav PPTA IPTA 
GWB’ 3.0 x 10715 1.45 x 10715 1.0 x 10715 1.7 x 10715 
Lentati et al. [90] Arzoumanian et al. [9] Shannon et al. [131] Verbiest et al. [153]? 

è - 1.5/yr 0.75/yr - 
BWM — Arzoumanian et al. [8] Wang et al. [154] — 
cwd 1.5 x 10714 3.0 x 10714 1.7 x 107! — 


Babak et al. [13]? Arzoumanian et al. [6] Zhu et al. [159] — 


^ Limits on the GWB are typically given as upper limits on the dimensionless strain amplitude at a 
frequency of 1/yr. 

b Verbiest et al. [153] note that the limit derived from the IPTA data set was only indicative and 
not rigorous as a full analysis was deferred to a future paper. 

€ BWMs can be quantified in many ways. In order to provide some comparative measure, this table 
presents the upper limit on the burst rate for bursts with normalised characteristic strain amplitude 
107", 

d Limits for CWs are given as upper bounds on họ, at a GW frequency of 10 nHz. 

° Babak et al. [13] don’t give a specific value; the number given is based on their Figure 3. 


Recent and Ongoing Improvements in PTA Sensitivity 


Much of the present interest in PTA experiments derives from the work by Jenet 
et al. [67] who predicted that a GWB should be detectable after a mere 5 years of 
timing on 20 pulsars, with a timing residual RMS of 100 ns — a level of precision 
that had only recently been demonstrated to be achievable in practice [141]. Sub- 
sequently, it became clear that the ideal PTA (20 pulsars, 5 years and 100 ns RMS) 
was unlikely to ever become a reality since the pulsar population is by nature highly 
inhomogeneous, implying a few pulsars would be likely to be timed at better pre- 
cision than 100ns, but most probably would not. Consequently, scaling relations 
were derived, first by Jenet et al. [68] and later by Siemens et al. [133], to allow 
fine-tuning of PTA experiments with a view to optimising sensitivity and shortening 
detection time scales. 

Siemens et al. [133] showed that the S/N of a GWB in a PTA data set scales with 
typical properties of the data set as follows: 


S/N x NCATE rg ts, (3) 
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where N is the number of pulsars, C is the cadence of the observations, A is the 
amplitude of the GWB, T is the length of the timing data set and o is the RMS of 
the timing residuals. While this analysis makes some basic simplifications in terms 
of data homogeneity, it does show clearly the very strong dependence of PTA sensi- 
tivity on the number of pulsars in the array. Consequently, a large number of pulsar 
surveys have been undertaken in recent years to increase the number of PTA-useable 
MSPs, as described below. The dependence on all other parameters is far less sig- 
nificant, except for the data length, which still adds considerably. In the low-S/N 
regime, Siemens et al. [133] show that a PTA’s sensitivity scales most strongly with 
the data length: S/N œ T'3/3, in close agreement with the earlier findings of Jenet 
et al. [67]. For long data sets the sensitivity is however affected by the level at which 
other long-period signals can be mitigated. Most significantly this refers to DM vari- 
ations which can be mitigated provided the observing set-up has been well chosen. 
Pulsar timing noise is equally important, but in the absence of predictive models or 
independent estimates of this noise source, the only possible approach is to limit 
its impact by post-facto modelling and subtraction. Finally, in the intermediate S/N 
regime (i.e. as we get closer to detections rather than mere limits), the sensitivity 
to a GWB is only weakly related to the timing precision of the MSPs, but for sin- 
gle sources of GWs the timing precision is still the dominant factor, consequently 
some efforts are being made to further improve the timing precision of MSPs by 
instrumental improvements and building of new telescopes. At the end of this sec- 
tion, a brief overview is given about various studies that quantify how all of these 
improvements are likely to affect the time to the first GW detection with PTAs — 
which mostly agree a detection within years to at most a decade is highly likely. 


Pulsar Surveys 


In order to increase the sample of MSPs that can be used in PTA experiments, a num- 
ber of pulsar surveys have been undertaken in recent years and are being planned 
for the near future. Specifically, two long-lasting surveys have been running on the 
Arecibo radio telescope for most of the past decade: the P-Alfa survey [37] and the 
Arecibo drift-scan survey [42]. In addition, the Effelsberg and Parkes radio tele- 
scopes (in Germany and Australia respectively) are continuing their all-sky partner 
surveys HTRU North and South [18, 73]. Parkes has simultaneously been equipped 
with cutting-edge processing technology as part of the SUPERB [72] survey, which 
aims to do real-time searches for pulsars and fast radio bursts. Data from the Green 
Bank Telescope (GBT) continues to be analysed in the Green Bank Northern Celes- 
tial Cap (GBNCC) survey [140]. At low frequencies, the LOFAR telescope [148] is 
finishing processing of the LOFAR Tied-Array All-Sky Survey [LOTAAS, 126] and 
at the GMRT the GMRT High-Resoluiton Southern Sky survey (GHRSS) is ongo- 
ing [23]. New telescopes are also getting up to speed on pulsar surveys, with the first 
successful discoveries published by the FAST telescope [158, 120], as part of the 
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Commensal Radio Astronomy FAST Survey (CRAFTS) and survey observations 
recently commenced for the MeerKAT “TRAPUM” survey [137]. 


ITSM Studies 


In order to increase sensitivity of pulsar-timing data sets to long-term signals like 
those expected from a GWB, it is of utmost importance to understand, mitigate and 
model any long-term signals that may be affecting the timing. Most such effects are 
deterministic effects contained in the timing model, but two more complex sources 
of red noise exist: timing noise and IISM noise. As discussed earlier, the origin of 
timing noise has not been unequivocally determined and consequently most models 
are rather ad-hoc power-law descriptions of uncorrelated timing signatures. ISM 
noise (or DM variations) is different, since it is the only known effect that causes a 
frequency dependence in timing residuals. 

The frequency dependence of ITSM noise implies that in principle it can be mea- 
sured and modelled independently from all the other timing-model parameters and 
correlated signals, because it can fully rely on the frequency resolution of the data. 
Specifically, three scenarios could be envisaged for measuring and correcting DM 
variations in pulsar-timing data [66]: 


e Multiple different observing bands: When more than one observing band is used, 
the frequency difference between the bands can be used as a lever arm that en- 
ables high-precision DM estimates. This idea has been implemented both with 
co-axial receivers and, more recently, with ultra-broadband receivers such as the 
UWL in Parkes [61] and similar observing systems at the Effelsberg and Green 
Bank radio telescopes, often not using actual instantaneous observations, but by 
determining an average DM over a range of dates [74]. In this scenario, care must 
be taken in the allocation of the observing time, since the DM sensitivity of the 
data will scale with the square of the observing wavelength, but the timing pre- 
cision may not, depending on the spectral index of the pulsar. Specifically, since 
on average pulsars have a spectrum that is flatter than v~? [65] and since the sky 
background noise is steeper than typical pulsar spectra [155], it is likely that the 
low-frequency band will still have superior sensitivity to DM, but have worse 
timing precision overall. For observations beyond 1.4GHz, however, the situa- 
tion often reverses, in that higher-frequency bands may have both less sensitivity 
to dispersion and lower timing precision [see, e.g. 83]. This implies a change in 
integration time depending on the observing frequency, may be in order. A more 
extensive analysis of post-correction ToA precision in this scenario, is given by 
Lee et al. [88]. 

e Low-frequency monitoring: As a possible way to mitigate the complexities of 
balancing DM and timing precision, one could attempt to monitor pulsars at low 
frequencies (generally at or below 400 MHz) and derive independent DM time 
series from those low-frequency data, to correct the higher-frequency data. This 
has the additional advantages that the DM modelling is now fully independent 
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from the higher-frequency timing; and at low frequencies DM corrections could 
often be measured within a single observation, which avoids correlations with 
effects like timing noise or other timing-model parameters. This approach has 
its own drawbacks because other IISM effects like scattering also become more 
pronounced at low frequencies and so the DM measurements may be biased or 
corrupted. Finally, in some cases the differences in Fresnel scale!’ at the top 
and bottom of the observing band make it possible that the actual space probed 
by electrons at different frequencies is slightly different — and hence the DM as 
well. This results in a frequency-dependent DM, which has been theoretically 
predicted [38] and observed [45], but the overall impact of this phenomenon on 
PTA sensitivity has so far not been accurately quantified. 

e High-frequency timing: Finally, with sufficiently sensitive telescopes, the main 
observing frequency could be moved to higher frequencies, where the IISM ef- 
fects are weaker. Pulsars also tend to be fainter at those frequencies, but in the 
case of highly sensitive telescopes, this may be a blessing since it implies jitter 
noise will be less significant, as more pulses will need to be averaged per obser- 
vation. However, with present telescopes, this approach may require too much 
observing time [83] — with more sensitive, future, telescopes it may become a 
realistic option. 


When it comes to effects of the ITSM, dispersion is only the peak of the iceberg. 
Things get a lot more complex when we consider multi-path propagation or scat- 
tering. This phenomenon arises when some of the photons get deviated from the 
straight line between the pulsar and Earth due to refraction; and later get refracted 
back into the line of sight. In its simplest form, for a thin scattering screen with 
density inhomogeneities that follow a Kolmogorov spectrum, this should cause a 
delay in photon arrival time which scales with the observing frequency to the fourth 
power: 

T scat CS SC 
where the scattering index @ is theoretically expected to be 4.4. In practice, Bhat 
et al. [21] have measured an average spectral index of the scattering strength of 
a = 3.9+0.2, slightly inconsistent with theory. More detailed studies at lower fre- 
quencies [51] have shown that in many cases the frequency scaling of the scattering 
was more consistent with a highly anisotropic scattering screen than with the typical 
Kolmogorov screen. 

Regardless of the frequency scaling, the primary observable effect of scattering 
in pulsar timing is that the pulse profile gets smeared out and gets a characteristic 
exponential tail. This worsens timing precision since it can wash out features and it 
may corrupt the DM measurement (although absolute DM measurements may not 
be necessary for pulsar-timing experiments anyway), but in principle this would not 
affect GW sensitivity as long as the effect is constant in time. 

The strength of scattering does change in time, though, as can most readily be 


seen by inspecting a “dynamic spectrum”: a plot of pulsed intensity as a function 


10 The Fresnel scale is a basic measure for the size of the Fresnel zone, which in turn is the region 
of space a signal can travel through between emitter and receiver. 
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of frequency and time. Such time-variable spectra often show changes in pulse in- 
tensity — a phenomenon known as scintillation. Diffractive scintillation (as shown 
in Figure 6) occurs when photons travelling from the pulsar to Earth meet refractive 
structures and get slight phase shifts due to location-dependent refractive indices. 
These phase-shifts cause constructive and destructive interference which are seen 
as bright and dark pathes in the dynamic spectrum. Since scattering and diffractive 
scintillation are effectively two different observables caused by the same turbulent 
and diffractive structure, it should not come as a surprise that they are related — in 
fact, they are inversely proportional [93]. An extensive analysis of PTA data car- 
ried out by Levin et al. [93] quantified the variations in diffractive scintillation and 
consequently estimated how variable scattering is in typical PTA observations. This 
analysis showed that at present levels of sensitivity, scattering variations are only 
rarely a real concern [with one exception studied in detail by 92], but in the next 
era of highly sensitive telescopes (notably MeerKAT, FAST and SKA), this will 
probably change. 

A further complication could arise from a higher-order effect where scattering 
and diffractive scintillation combine in what are known as “scintillation ares” [139]. 
This primarily occurs in highly anisotropic media, when the majority of the radi- 
ation comes from the pulsar (essentially a point source), but significant amounts 
of energy come from other, typically straight and narrow, structures on the sky. It 
causes ripples across the dynamic spectrum, which are more easily noticed in the 2D 
Fourier transform of the dynamic spectrum — this is also called the “secondary spec- 
trum” (see right-hand side of Figure 6). The initial discovery of secondary spectra 
is still relatively recent and since these tend to be faint features which require high 
sensitivity and high resolution (i.e. large data rates), their study has only developed 
slowly. However, early studies by Hemberger and Stinebring [58] already showed 
how these phenomena can impact timing significantly. Turning this around, Reardon 
et al. [121] showed how scintillation arcs can actually be useful for timing, as they 
can provide independent constraints on timing-model parameters. The study of how 
to use scintillation arcs, or what effects they really do have on PTA experiments, has 
only just begun, so at present their impact is not fully clear yet. 

A final concerning occurrence that the ITSM may create, are frequency-dependent 
DMs (also often — and confusingly — named “chromatic DMs”). The principle be- 
hind frequency-dependent DMs is as follows: in wide-band observations, photons 
with a wide range in wavelengths are observed. These photons did not all travel 
through the same space — in fact, since the Fresnel scale is frequency-dependent, 
there is a bias that causes lower-frequency photons to be able to travel through a 
wider region of space than higher-frequency photons. If the variations in electron 
density are sufficiently extreme — or if measurement precision is sufficiently high — 
this would imply that the high-frequency photons may sample a different electron 
distribution than the low-frequency photons. Consequently, the DMs measured at 
the top and bottom of the observing bands may be different because they refer to 
different parts of space. 

While the concept of frequency-dependent DMs had been known for much 
longer, the theoretical description was first laid out by Cordes et al. [38]. Detec- 
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Fig. 6 Dynamic and secondary spectrum of PSR J0826+2637 (PSR B0823+26). Left: the dynamic 
spectrum of a half-hour observation on PSR J08216+2637 with the LOFAR core. Shown is the 
pulsed intensity on an arbitrary colour scale (units are uncalibrated) for a segment of S-MHz band- 
width, centred on 147.5 MHz. The bright patches are called scintles and are caused by diffractive 
scintillation. Less clear is the higher-frequency corrugations that run diagonally across this dy- 
namic spectrum and which are caused by a combination of diffractive and refractive scintillation. 
Right: this figure shows the secondary spectrum of the left-hand plot, i.e. the 2D Fourier transform, 
whereby the power levels are shown on a logarithmic scale. A highly asymmetric arc is visible, 
extending out to fractions of a millisecond at negative fringe frequencies (also called “Doppler 
rates”) but only out to about 0.1 ms at positive fringe frequencies. Figure courtesy of Ziwei Wu. 


tion of such a phenomenon is naturally complex given the many other frequency- 
dependent effects described earlier, but by looking at the time-difference of DMs 
measured at opposite parts of a low-frequency observation, Donner et al. [45] suc- 
ceeded in making the first clear detection of such chromatic effects. These initial re- 
sults showed rather a more complex picture than the theory had predicted, clarifying 
that further research into frequency-dependent DMs is required before a conclusive 
understanding of their possible impact on PTAs can be drawn. Given the continuous 
coverage over extremely wide frqeuency ranges of new telescopes like the uGMRT, 
ngVLA, MeerKAT and the SKA, a much clearer understanding is bound to arise 
within the next decade. 


Sensitivity Predictions 


As mentioned earlier, as PTAs edge closer to a GW detection, the number of pulsars 
in the array is of key importance for the PTA’s sensitivity. Since predicting pulsar 
discoveries is infamously hard, it is equally hard to make accurate predictions of 
PTA sensitivity given future pulsar discoveries, yet several papers have demontrated 
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the validity of the adage “More pulsars is more sensitivity”. Most recently, Kelley 
et al. [75] showed that — assuming ongoing regular discoveries of MSPs that can be 
timed at high precision — all PTAs could hope to detect GWs within a few years; 
and all were virtually guaranteed a detection within a decade. 

Based on the most up-to-date predictions for a GWB in the PTA band and real- 
istic numbers from existing PTA experiments, Rosado et al. [124] showed that the 
most likely scenario would be that a GWB would be detected in about one to two 
decades time. A detection of CWs was also not unrealistic, but would probably take 
somewhat longer. Rosado et al. [124] did not investigate the impact of increasing 
the number of MSPs in the array, but did evaluate the impact of more sensitive sys- 
tems — the SKA in particular and found that with the full SKA, a GWB should be 
detectable within a few years; and CWs within about a decade. 

The idea that a highly sensitive telescope could detect GWs within a few years, 
even with existing pulsars, was also demonstrated by Lee [85], who drew essentially 
the same conclusion for the CPTA with its unprecedentedly sensitive set of tele- 
scopes. Also Lazarus et al. [83] investigated the impact of improved sensitivity on 
PTA detection time scales, this time in the context of improved receiver and record- 
ing systems. Specifically, they estimated that the new recording system at Effels- 
berg would improve the telescope’s sensitivity to a GWB by a factor of up to three 
compared to the status quo, in only four years time. They furthermore continued the 
work by Lee et al. [88] to demonstrate how wide-band and low-frequency observing 
systems might aid the detection of GWs by efficiently correcting DM variations and 
thereby keeping the timing RMS low. Verbiest et al. [153] approached the sensitivity 
improvements in a very complementary way, showing that global collaboration and 
sharing of data would lead to approximately a factor of two improvement in GW 
sensitivity. 

It is also possible that the coming years continue to bring non-detections; this 
could occur if we encounter an unexpected instrumental noise floor, such as intrin- 
sic pulsar noise or a high level of ephemeris uncertainty. However, this scenario is 
unlikely to be an issue. Pulsar noise can be overcome by targeted noise modelling 
[as in 55], or by simply adding more pulsars to a PTA, which will beat down noise 
that is uncorrelated between different pulsars, thus still raising our sensitivity to the 
correlated GW signals. Regarding uncertainties in the Solar-System ephemerides, 
as previously noted on page 15, ephemerides uncertainties are correlated. The sig- 
nal, however, is dipolar thus we can attempt to measure and remove it, even if 
there is some leakage between dipolar and quadrupolar signals [145]. In addition, it 
has been demonstrated that PTAs are already breaching the accuracy of published 
ephemerides, and techniques have been developed to overcome such uncertainties 
[147]. 

Regardless, current upper limits on the GWB are already impacting our under- 
standing of the evolution of galaxies and their supermassive black hole residents 
[131, 134, 9, 34], in addition to placing novel constraints on cosmic strings [9, 156] 
and exotic forms of dark matter [e.g. 35]. For the interested reader, the wealth of ac- 
cessible science with GW limits and detections is the subject of another large review 
[30]. 
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For the purposes of this review, it suffices to say that if our limits continue to 
improve around an order of magnitude beyond their current point, it would be as- 
trophysically surprising. This is because even the most pessimistic simulations of 
supermassive binary black hole evolution (where no galaxy merger results directly 
in a binary supermassive black hole due to inefficient inspirals) still result in signals 
detectable at the hy; > 107 16 level [25]. See also the “Massive Black Hole Mergers” 
chapter of this edition by E. Barausse and A. Lapi. 

In summary, there are a variety of ways in which PTAs can further gain sensi- 
tivity and all of these ways are being explored. Essentially all predictions conclude 
a detection is likely within the next few years or at most within the next decade — 
regardless of which particular improvement is being studied. This is not surprising 
given that our most reliable predictions on the strength of the GWB are right up 
against our most constraining limit [131, 9], which strongly implies a detection is 
bound to be imminent. 


Summary 


The extreme properties of the spinning neutron stars called pulsars enable a unique 
experiment to detect GWs in a spectral range that is highly complementary with 
other mature GW projects like LIGO and LISA. PTAs are expected to make the 
first detection of nHz GWs within the next few years and in doing so, will allow 
new and unprecedented constraints on galaxy formation and evolution scenarios. 
The way forward is long and hard, however, as data sets are complex and highly 
heterogeneous and a variety of noise sources need to be dealt with. Instrumental 
upgrades and extension and improvement of observing schedules and source lists 
are underway to further enhance sensitivity — a process that will culminate in the 
ultimate PTA to be ran on the telescope of the future: the SKA. With the added 
collecting area and the larger number of pulsars that the SKA will be able to time 
at high precision, GW astronomy in the nHz range can be expected to properly take 
off. 


Cross-References 
Barausse & Lapi, “Massive Black Hole Mergers” 
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